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Abstract. This work investigates the possibility of particle-based algorithms for the Navier-Stokes equations and higher order 
continuum approximations of the Boltzmann equation; such algorithms would generalize the well-known Pullin scheme for 
the Euler equations. One such method is proposed in the context of a discrete velocity model of the Boltzmann equation. 
Preliminary results on shock structure are consistent with the expectation that the shock should be much broader than the 
near discontinuity predicted by the Pullin scheme, yet narrower than the prediction of the Boltzmann equation. We discuss the 
extension of this essentially deterministic method to a stochastic particle method that, like DSMC, samples the distribution 
function rather than resolving it completely. 
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INTRODUCTION 

In computing flows containing both rarefied and continuum regions [1], information transfer between regions is 
simplified if both are computed by a comparable algorithm. Hybrid particle -based Euler-DSMC has been demonstrated 
[2] using the Pullin scheme for the Euler region. However, a smooth transition between the solution schemes is best 
effected if an overlap region exists in which both methods are valid [3]; since the Euler equations are the zero Knudsen 
number limit of the Boltzmann equation, the existence of such an overlap region for hybrid Euler-DSMC in general 
might seem problematic. This observation motivates developing particle methods for Navier-Stokes, and possibly for 
higher order continuum approximations as well [4]. This paper has the restricted goal of obtaining Navier-Stokes 
dynamics using a discrete velocity model (DVM) [5] of the Boltzmann equation. We outline how the method might 
be applied to a stochastic particle implementation closer in spirit to DSMC, in which the distribution is sampled by a 
small number of particles rather than resolved completely, and briefly consider some aspects of the numerical diffusion 
in such particle methods for continuum flows. 


DISCRETE VELOCITY MODEL 

As in the formulation of Gatignol [5], we consider a finite discrete velocity set {c,}, which has the form {c,} = 
U{(±«i, ±« 2 ,±« 3 )} where the n,- are integers; invariance under the 24 symmetries of the cube is assumed. All discrete 
velocities such that e(n\,n 2 ,n 3 ) = n\ + n\ +n \ < £ max are chosen, so that the velocities form a ‘lumpy sphere.’ In a 
model with £ max = 4, there are 33 discrete velocities; in a model with £ max = 9, there are 123. But unlike the models 
discussed in [5], our DVM is implemented on a lattice consistent with the discrete velocity space, so that c,-Af is an 
integer multiple of the lattice spacing Ax. Therefore, ‘molecular’ motion always carries particles from one lattice site 
to another. This imposes a key property of the Boltzmann equation, that collisions occur at a point, and eliminates 
numerical diffusion. Let N denote the number of particles with velocity c so that the distribution function /(x,c) 
takes the discrete form }_N l (\)8(c — c then the discrete Boltzmann equation is taken to be 

Ni{x + aAt,t+At) -Ni(x,t)=At Y, £lij Pq {-Ni{x)Nj(x) +N p (x)N q {x)} (1) 
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where all collisions consistent with momentum and energy conservation are performed. This number is very large: 
for 33 velocities, 3,288 collisions are possible, and for 123 velocities, 96,027 collisions are possible. The factor Q. is 
defined with a discrete version of hard-sphere collisions, but no effort is made to model realistic collision dynamics. 
It might be appropriate to call this model a ‘lattice discrete velocity model.’ While it does represent a possible particle 
dynamics and can therefore mimic many properties of the Boltzmann equation, the finite velocity space imposes some 
special requirements which will be noted subsequently. 



In what follows, it will be convenient to use notation appropriate for the continuous case when the extension 
to the discrete case is trivial, and to describe changes required in the discrete setting only when necessary. Thus, 
given the distribution function /, consider f MB , the equilibrium distribution with the same hydrodynamic moments 

as /. We have the usual conserved hydrodynamic moments p = J fdc = J f MB dc, pu = j \fd\ = J \f MB d\, 
MB c 2 dc, here, v is the particle velocity, u is the hydrodynamic velocity, and c = v u is the peculiar 
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velocity. A property of this DVM is that the equilibrium is formally Maxwellian; however, because of the finite velocity 
space, the parameters in the Maxwellian need not coincide with the conserved moments. This fact forces us to evaluate 
f MB by relaxing / to equilibrium by multiple collision steps, even if this process is time-consuming. 


DISCRETE FORMULATION OF THE CHAPMAN-ENSKOG EXPANSION 


The discrete Boltzmann equation, Eq. (1), is solved by the usual advection-collision operator splitting. Beginning 
with an approximate distribution function /„, let / 1 = /„ + (Af)c • (the advection step), and f n+ \ = f 1 + 
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c • +/[/o,/o] and we are done; otherwise, these steps are 
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(A t)J[f n+ i,f n+ i\ where J denotes the collision integral defined in Eq. (1) (the collision step). If f n+ \ — f n = 0(At 2 ) 

then fn+\ ~ fn = fn + (Af)/[/ n+ 1 ,f n+ 1 ] - /o = (At) 
repeated. 

The Pullin scheme is a simple modification in which f n+ \ = f M, \ . Thus, the collision step is replaced by instanta- 

n-\- ^ 

neous relaxation to a Maxwellian. But as noted before, in the DVM used here, f n+ \ is evaluated by repeated collisions: 
in the problem discussed later, four iterations led to an equilibrium distribution, but no general statement can be made 
about the number of iterations necessary. As before, we iterate until f n+ \ = /„ to order At 2 . If this actually happens, 
then f n+ 1 satisfies the steady Boltzmann equation; but since f n+ \ is Maxwellian, the hydrodynamic moments satisfy 
the Euler equations. 

As stated above, our problem is to go beyond Euler, to higher order continuum approximations. We begin by consid- 
ering the Grad 13-moment expansion [6] / = /-°* +/ ^ ^ where f^ = f MB . (Recall that/^, the projection along 

c, vanishes.) Grad’s results f ^ = jf^a^ : H^ 2) and / ( ' 3 ’ = -p)/^ 0l b^ 3) : H (3) both apply in the DVM (we distinguish 

tensors and vectors by font and use index notation only where necessary to avoid ambiguity), where H [2) (c) = 3 — cc — I 


(c is the peculiar velocity defined above, and I is the 


identity), and H (3 )(c) = 
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c. However, as 


noted previously, the finite velocity space imposes some modifications: one is that in the second-order Hermite poly- 
nomial, I = - J t/c/ ,: 0 : cc: the continuous result Ijj = S ;/ is not exactly true in the DVM. One reason is that the discrete 

velocity space is at most invariant under the finite group of rotations of a cube, not the continuous group of rotations 
of space. Therefore, invariant tensors of second rank can exist other than the unit tensor; some of the implications are 
developed in detail in the context of lattice Boltzmann models in [7]. Another modification imposed by the discrete 

velocity space comes about because the definition = — / dc /H' 2 ' from the continuous Grad expansion cannot be 
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used. The reason is that the crucial property J t/c{l,c, c 2}y(p) = o for p > 1 which insures that the hydrodynamic 

moments are carried by the lowest order approximation f(°\ is not satified in the discrete setting. This fact forces a 
modification for the DVM, namely that / ^ : H^ 2 ) where 


b (2) = a (2)_ a (2) :H (2) 


H(2) : H(2) 
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Note however, that the magnitude of such corrections decreases if the number of discrete velocities is increased. The 
third order term pertaining to the heat flux vector is treated similarly, but for brevity, explicit expressions are not given 
here. 

We propose to implement the 13-moment approximation in this DVM following the Pullin scheme, by defining 
/„ + i = f'^l , + /^i + f^i at each collision step. Thus, the distribution f n+ i resulting from the advection step is 



first relaxed to equilibrium, to give / 0 , as in the Pullin scheme. Next, this distribution function is projected onto 

( 2 ) ( 3 ) 

the relevant Hermite polynomials of the Grad expansion to give / , and f , . Because the Grad expansion is purely 

local, these steps are all straightforward. If iteration of these steps converges, the result will be an approximate solution 
of the Boltzmann equation by a Grad 13-moment distribution. It is perhaps more precise to say that if the scheme 
converges, it solves the Boltzmann equation to the order of the Grad 13-moment equations: because the distribution 
function evolves freely between collision steps, its evolution may not be determined entirely by the evolution of its 
moments (compare [4]). 

This procedure was tested by applying it to the one-dimensional problem of particles with a constant hydrodynamic 
velocity impinging on a specular wall, here modeled by the ‘bounce-back’ boundary condition of the lattice Boltzmann 
equation. The particles are confined to a finite interval, and therefore separate from the left hand boundary, forming 
a rarefaction wave, and are reflected from the right hand boundary, forming a shock wave. Figure 1 compares the 
unmodified DVM to the Pullin and Grad schemes just described. The most important feature is that the reflected shock 
moves at the same speed in all three cases; the use of Eq. (2) is crucial to this property, otherwise, /'■ 2 1 incorrectly 
contributes to the total energy and makes the shock move faster. At this point, we only assert qualitative agreement 
with the expectation that the Grad scheme shock is less steep than the shock predicted by the Pullin scheme, but steeper 
than the very broad shock predicted by the unmodified model. 
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In principle, the transition from Grad to Navier-Stokes only requires introducing the approximations a u = 
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and a 1 ; 1 ’ 1 = 


the definitions of /( 2 ) and but this is disadvantageous 
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because it requires arbitrary choices in the evaluation of partial derivatives of the hydrodynamic variables; it is prefer- 
able to express everything in terms of particle operations. To this end, we recall that [8] 
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where & = — + v • — is the advection operator, and C ( ° is the trace-free part of CjCj, and similarly 


But we obviously have [8] 
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where a and p depend on the collision model. These functions could be evaluated for any given model, including the 
DVM; in the continuous case of Maxwell molecules, we have, in the notation of Waldmann [9] 
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where the subscript M denotes the linearized collision operator for Maxwell molecules, and ft>n = (p/ p )(d\ \ with 
a»n =2/3 and ©20 = (p/jtt )© 20 with ©20 = 1. It follows from Eqs. (3) and (4) that for Maxwell molecules. 


( 2 ) = P ©11 g (Q) _ / P 

p©n-©20 ' \P J ffln-©20 


1 




(5) 


with a similar result for ; for any molecular model, the corresponding result is found by solving a trivial system 
of linear equations. We therefore have the required result: an expression for a Chapman-Enskog distribution in terms 
of advection and collision operations alone. Replacing f' 2 - in the discrete Grad scheme by Eq. (5) and f' 22 ' 1 by the 
easily derived analog, leads to a method which, if it converges, converges to a solution of the Boltzmann equation to 
Navier-Stokes order. 



PARTICLE IMPLEMENTATION AND NUMERICAL DIFFUSION 


The Grad and Navier-Stokes particle schemes just described are deterministic. An implementation closer in spirit 
to DSMC is obtained by locating a finite number of particles at each lattice site, advecting them by their discrete 
velocities between lattice sites, and colliding them by stochastic sampling of the possible collision outcomes. As in 
DSMC, this procedure evolves a stochastically sampled distribution function, rather than the distribution function 
itself. The moments of the sample can be used to evaluate an approximate Maxwellian f MB at every point as in the 
original Pullin scheme; since it is expressed in terms of the advection and collision operations, it should be possible 
to evalute an approximate Chapman-Enskog distribution as well and thus to construct a solution of the Navier-Stokes 
equations entirely by a particle method. 

This type of continuum flow solution will be useful only if it can be applied to regions larger than the mean free path, 
but this raises the problem of numerical diffusion. The problem can be illustrated in the DVM by applying collisions 
to neighboring lattice sites, so that the term 

Yj &ij Pq {-Ni{x + Ax)Nj{x) +N p (x + Ax)N q (x) -Ni{x)Nj{x + Ax) + N p {x)N q {x + Ax)} (6) 
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and the same term with —Ax replacing Ax, are added to the collision term in Eq. (1). The results are illustrated in 
Figure 2. The graph on the left shows the shock region in the problem treated in Figure 1 computed by the unmodified 
DVM. This result was compared to the result of following each collision step with collisions between adjacent lattice 
sites following Eq. (6). 

The perhaps surprising result is that the shock front steepens, but this was due to the effect of additional collisions. 
A better comparison is with the result of performing two collision steps at each lattice site for each advection step. This 
comparison shows that collisions with adjacent sites produces a small but noticeable broadening of the shock front, 
which can be considered the result of diffusion due to collisional (in additional to advective) coupling between different 
sites. The graph on the right makes the same comparison for the Pullin scheme: the Pullin scheme is compared to a 
scheme with collisional coupling between adjacent sites with the same number (8 in this case) of collision steps per 
advection step. In this case, the diffusion effect is extremely small, perhaps even surprisingly small, but it is definitely 
nonvanishing. It is hoped that the DVM can provide a useful way to test ideas for mitigating this numerical diffusion. 

A possible model for collisional coupling between distribution functions at different points is the diffusive Boltz- 
mann equation 




where v is an empirical diffusivity. Under this model, the hydrodynamic equations including the continuity equation, 
all acquire diffusion terms. But in the Chapman-Enskog expansion, the added terms appear in a different way from the 
standard ones, because in effect, one factor of c is replaced by vd jdx. The result is coupling between the velocity and 
temperature diffusion. 


CONCLUDING REMARKS 

It appears to be possible to generalize the Pullin scheme so that it computes a solution of the Navier-Stokes equations. 
This is accomplished by adding suitable higher order moments to the Maxwellian after each advection step and 
modifying the moments so that they are proportional to gradients of the hydrodynamic variables. The important fact is 
that these gradients can be computed in terms of the particle operations of advection and collision; it is not necessary to 
introduce explicit numerical differentiation. It remains to be shown that these steps do not require complete resolution 
of the distribution function, but that they remain consistent with stochastic sampling of the distribution function, as in 
particle methods like DSMC. The effective diffusion caused by particle collisions over distances larger than the mean 
free path remains an important limitation. 
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FIGURE 1. Particles are moving at constant hydrodynamic velocity toward the specular right-hand wall. LEFT: the complete 
density field, comparing the ‘Boltzmann’ (red solid line), ‘Euler’ (dashed green line), and ‘Grad’ (blue dot-dash line) solutions as 
described in the text. RIGHT: Same key as on the left, the reflected ‘shock’ region alone is shown at two different times, to verify 
agreement of the front velocity in all three models. 
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FIGURE 2. LEFT: shock region, unmodified DVM (red solid line), collisions with nearest neighbor lattice site added (dashed 
green line), and DVM with two collision steps per advection step (blue dot-dash line) RIGHT: shock region, Pullin scheme (blue 
dot-dash line) and Pullin scheme with collisions with nearest neighbor lattice site added (dashed green line) - note that the x scale 
is much finer than on the other graphs 
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